Bubble merging in breathing DNA as a vicious walker problem in opposite potentials 
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We investigate the coalescence of two DNA-bubbles initially located at weak domains and sepa- 
rated by a more stable barrier region in a designed construct of double-stranded DNA. In a contin- 
uum Fokker-Planck approach, the characteristic time for bubble coalescence and the corresponding 
distribution are derived, as well as the distribution of coalescence positions along the barrier. Below 
the melting temperature, we find a Kramers-type barrier crossing behavior, while at high tempera- 
tures, the bubble corners perform drift-diffusion towards coalescence. In the calculations, we map 
the bubble dynamics on the problem of two vicious walkers in opposite potentials. We also present 
a discrete master equation approach to the bubble coalescence problem. Numerical evaluation and 
stochastic simulation of the master equation show excellent agreement with the results from the 
continuum approach. Given that the coalesced state is thermodynamically stabilized against a state 
where only one or a few base pairs of the barrier region are re-established, it appears likely that this 
type of setup could be useful for the quantitative investigation of thermodynamic DNA stability 
data as well as the rate constants involved in the unzipping and zipping dynamics of DNA, in single 
molecule fiuorescence experiments. 

PACS numbers: 02.50.Ey, 05.40.Fb, 82.20.-Uv, 82.37.-j, 87.14.G- 



I. INTRODUCTION 



Within a broad range of salt and temperature condi- 
tions, the Watson-Crick double-helix [l| is the equilib- 
rium structure of DNA. This thermodynamic stability is 
effected by hydrogen-bonding between paired bases and 
by base-stacking between nearest neighbor pairs of base 
pairs [H, i, i, [1 i, 0, 0, Hi- By an increase of the tem- 
perature or by variation of the pH-value (titration with 
acid or alkali) double-stranded DNA progressively dena- 
tures, yielding regions of single-stranded DNA, until the 
double-strand is fully molten. This is the helix-coil tran- 
sition ■ The melting temperature Tm is defined as 
the temperature at which half of the DNA molecule has 
undergone denaturation 0, i, [m, m. Typically, the de- 
naturation starts in regions rich in the weaker Adenine- 
Thymine base-pairs, and subsequently moves to zones of 
increasing Guanine-Cytosine content. The occurrence of 
zones of different stability within the genome was shown 
to be relevant when separating coding from non-coding 
regions [il,[l3. 
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However, already at room temperature thermal fluc- 
tuations cause rare opening events of small intermittent 
denaturation zones in the double- helix [l^ . These DNA 
bubbles consist of flexible single-stranded DNA, and their 
size fluctuates by step-wise zipping and unzipping of the 
base pairs (bps) at the zipper forks, where the bubble 
connects to the intact double-strand. Initiation of a bub- 
ble in a stretch of intact double-strand requires the cross- 
ing of a free energy barrier AGbubbic of some 8 to 12 ksT 
at physiological temperature, corresponding to a Boltz- 
mann factor, often referred to as the cooperativity fac- 
tor, (To = exp(-AGbubbio/fcsr) - 10^^ . . . 10-3. Once 
formed below the melting temperature Tm, a bubble will 
eventually zip close. Above Tm, a bubble will preferen- 
tially stay open and, if unconstrained, grow in size until it 
merges with other denaturation bubbles, eventually lead- 
ing to full denaturation of the double-helix. Constraints 
against such full unzipping could, for instance, be the 
build-up of twist in smaller DNA-rings or the chemical 
connection of the two strands by short bulge-loops, com- 
pare Ref. [l6j . 

Biologically, the physical conformations of DNA 
molecules are considered of increasing relevance for its 
function, see, for instance, the review and refer- 
ences therein. In particular, the existence of intermittent 
(though infrequent) bubble domains is important, as the 
opening up of the Watson-Crick base pairs by breaking 
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of the hydrogen bonds between complementary bases dis- 
rupts the helical stack. The flipping out of the ordered 
stack of the unpaired bases allows the binding of specific 
chemicals or proteins, that otherwise would not be able 
to access the reactive sites of the bases 0, [l^ . In 
fact, there exists a competition of time scales between 
the opening/closing dynamics of DNA-bubbles and the 
binding kinetics of selectively single-stranded DNA bind- 
ing proteins [H, [l^, [10, Hlj . That the chemical potential 
of the single-stranded binding proteins does not lead to 
full denaturation of the DNA is due to the slow binding 
of the pro teins when compared to the bubble dynamics 
EliiSIIli . It is also beheved that DNA-breathingas- 
sists in the transcription initiation process pi. l23l. [23 . T25j . 
The quantitative knowledge of the denaturation dynam- 
ics as well as energetics is imperative to a better under- 
standing of genomic biochemical processes. 

DNA-breathing has been modeled extensively in terms 
of the Peyrard-Bishop model, that is based on a set of 
Langevin equations for the base-base distance in a base 
pair; the effective attraction between the bases is repre- 
sented by model potentials [H, [13, [H, [H, [s^l. Alter- 
natively, DNA-breathing can be considered as a random 
walk process in the free energy landscape of the Poland- 
Scheraga model of DNA denaturation, as the number of 
broken base pairs turns out to be the slow variable of 
the process |16| . In continuum form, this approach to 
DNA bubble dynamics has been described in terms of 
a Fokker-Planck equation [3l|, [13, [11]. A discrete de- 
scription, in which the coordinate of the random walker 
corresponds to a specific base pair, was suggested in 
Refs. |2ll [13, [H, [l^l, and the corresponding stochastic 
simulations analysis of DNA breathing was introduced 
in Ref. [13] ■ The influence of a random energy land- 
scape on bubble localization and dynamics was studied 
in Refs. [32, [13, [13], while a framework to include an 
arbitrary given sequence of base pairs was developed in 
Refs. [13, [H, [13]- Endowed with the sensitivity of their 
dynamics, DNA constructs were proposed as nanosensors 
[40l . [4]| . We note that the formulation in terms of the 
gradient of the Poland-Scheraga free energy allows one 
to explicitly introduce all necessary independent stack- 
ing parameters based on the study [3|, see also the dis- 
cussion in Ref. [13]. Measuring the dynamics of DNA 
bubbles also provides information on the magnitude of 
the critical exponent c representing the entropy loss fac- 
tor of a closed polymer loop, deciding on the order of 
the denaturation transition 0, \M, [H, [M [13, [H, [13 > by 
influencing the temporal survival probability of bubbles 

MM- 

The multi-state nature of DNA breathing can be mon- 
itored in real time on the single DNA level by fluores- 
cence correlation techniques [l6| . It has been shown in a 
quantitative analysis that the experimentally accessible 
autocorrelation function is sensitive to the stacking pa- 
rameters of DNA [13, [13]- However, it has not been fully 
appreciated to what extent the fluorophore and quencher 
molecules, that are attached to the DNA construct in 



the experiments reported in Refs. [13, [H, [13] , influence 
the stability of DNA. Moreover, the zipping rates mea- 
sured in the single molecule fluorescence setup differ from 
those determined in NMR experiments [13, [13] • We here 
propose and study a complementary setup for the sin- 
gle molecule fluorescence investigation of DNA breath- 
ing, as shown in Fig. [TJ In this setup, a short stretch of 
DNA, clamped at both ends, is designed such that two 
soft zones consisting of weaker AT-bps are separated by 
a more stable barrier region rich in GC bps. For simplic- 
ity, we assume that both soft zones and barrier are ho- 
mopolymers with a bp-dissociation free energy AGs and 
AGb, respectively, and, in accordance with the experi- 
mental findings of reference [l3], we neglect secondary 
structure formation in the barrier zone. At temperatures 
higher than the melting temperature T, of the soft zones, 
but still lower than the melting temperature Tt of the 
barrier region such that two open bubbles are being pro- 
moted, thermal fluctuations will gradually dissociate the 
barrier, until the two bubbles coalesce. Note that the 
melting temperature at 100 mM salt conditions differs 
by approximately 50 degrees between mixed (AT/TA)„ 
and (GC/CG)„ homopolymers respectively [3, [I]- This 
should provide a large enough temperature interval be- 
tween hard and soft zones to perform this type of ex- 
periment. In the following we use realistic values for the 
simulations. Once coalesced, the free energy correspond- 
ing to one cooperativity factor uo ~ 10^^ . . . 10^"^ is re- 
leased, stabilizing the coalesced bubble against reclosure 
of the barrier. This fact should allow for a meaningful 
measurement of the coalescence time in experiment, and 
therefore provide a new and sensitive method to measure 
DNA stability data and base pair zipping rates. We also 
study the case when the system is prepared as above and 
then T suddenly increased such that T > Tf, > Ts so 
that the system is driven towards coalescence. In both 
cases the two boundaries between bubbles and barrier 
perform a (biased) random walk in opposite free energy 
potentials. 

In fact, the study of the bubble coalescence is of in- 
terest in its own right, as we map the random walk of 
the two zipper forks separating double-stranded barrier 
base-pairs from already denatured single-stranded bub- 
ble domains onto a new case of the vicious walker prob- 
lem. Namely, we deal with two vicious walkers in lin- 
ear but opposite potentials. The viciousness condition 
corresponds to the fact that when the two zipping forks 
meet, the bubbles coalesce, and the dynamics is stopped. 
While the problem of a general number of (otherwise non- 
interacting) vicious walkers in free space was solved long 
time ago [50[ and has been only relatively recently gen- 
eralized to the case of motion in a common potential 
[Slj . even two walkers in different potentials cannot be 
addressed analogously by the straightforward antisym- 
metrization procedure. To solve our problem in the con- 
tinuum limit, we use a trick of introducing individual 
symmetry transformations for each walker which trans- 
form the respective Fokker-Planck operators to the same 
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Figure 1: Schematic picture of the bubble coalescence setup in a designed DNA construct. It is clamped at both ends and 
consists of two outer soft zones (thin red lines) of lengths Nl , Nr bps with melting temperature Ts and a stronger A'^-bps-long 
barrier zone (thick blue lines) with Tt > Ts- a) All bps closed {T < Ts < Tt,). b) Soft zones open by raising the temperature 
above Ts. bi-ba) Successive opening of the barrier driven mainly by fluctuations (T < Tb) or drift (T > Tt) until coalescence. 
The discrete coordinates X,Y = —Nl, . . . ,N + Nr are defined as the positions of the interfaces between the closed and broken 
bps. 



Hermitian form. In the transformed frame the problem is 
solved by the standard procedure constructing the joint 
probability density from the antisymmctrized product of 
the single- walker probability densities. This solution en- 
ables us to effectively reduce the numerical efforts needed 
for the evaluation of the joint probability density, mean 
coalescence time, spatial probability density of coales- 
cence position, etc. by one dimension. Moreover, some of 
the quantities of interest can be obtained this way fully 
analytically including all the characteristics of the high- 
barrier case which are very hard to reliably determine 
numerically. 

The paper is organized as follows. After introducing 
the discrete model in Sec. and its continuum limit 
in Sec. IIIII we solve the latter in Sec. |TV] by using the 
symmetry of the problem. The main quantities of inter- 
est, namely, the coalescence time and its distribution, as 
well as the distribution of the coalescence position are 
obtained in Sec. |Vl In the following Sec. IVII wc intro- 
duce a direct solution of the discrete problem via the 
complete master equation and a stochastic simulations 
scheme (Gillespie) and compare these results to those of 
the continuum approximation in Sec. lVIll In Sec. lVIlIl we 
address the connection between our models and results 
and real biologically relevant data. In the last Sec. lIXI we 
state our conclusions. Appendix |^ presents a detailed 
calculation of the auxiliary single-walker density. In Ap- 
pendix |B] we explain the direct numerical solution to the 



full master equation of Sec. IVII while in Appendix [Cl the 
Gillespie stochastic simulation scheme for the same mas- 
ter equation is briefly summarized. 



II. DEFINING THE ZIPPING AND UNZIPPING 
RATES 

In this section we define the transition rates for open- 
ing or closing a base pair, which are determined by two 
effects: the energy landscape stemming from the thermo- 
dynamical partition factors and the thermal fluctuations. 

As illustrated in Fig. [Tl we consider the case when the 
two soft bubble zones of the DNA construct are prefer- 
entially open, while the central barrier region is initially 
completely closed. We assume that the two soft zones 
are homopolymers with identical melting temperature Tg, 
and that the barrier region is a homopolymcr with melt- 
ing temperature Tf, > Tg. In the following, we neglect 
secondary structure formation in the bubbles, consistent 
with experimental observations in relatively short bubble 
domains p^ . The barrier region of initially closed base- 
pairs between the zipper forks will also be referred to as 
the clamp. 

Each of the two DNA-bubbles is characterized by a 
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partition factor of the form 

X = -Nl 

N+Nr 

= + + n <n (2) 

where X,Y ^ [—Nl, N + Nji] denote the positions of the 
left and right zipper fork, respectively. In Eqs. ([T]) and 
^ the quantity u(X), which takes the values Ms(b), is the 
Boltzmann factor for breaking a base-pair in the soft zone 
(barrier domain). Ugif,) = exp (/3AGs(b)) ? corresponding 
to the free energy AG's(b) for breaking a base-pair; more- 
over, P = l/[kBT]. We define u{~Nl) = uiN + Nn) = 1. 

Note once more the prefactor {bubble length + 1)^'^ dis- 
playing the inherent long-range character of the Poland- 
Scheraga free energy model. It measures the reduction of 
the degrees of freedom of a loop configuration, as char- 
acterized by the critical exponent c [3, EES [Hi. 
For the long-time behavior in larger, single bubbles the 
influence of c on the distribution of bubble lifetimes is 
considered in Refs. [s^ . [s^ in a continuum approach. 

Finally, (Tq is the cooperativity factor corresponding 
to the free energy barrier for breaking the first base-pair 
in a stretch of intact double-strand. Loosely speaking, 
it corresponds to the disruption of two stacking interac- 
tions in the DNA, while the single open base-pair's en- 
tropy gain cannot balance the required enthalpy. This 
contrasts the opening of further base-pairs, for which the 
entropy gain almost balances the enthalpy cost. The co- 
operativity factor (Jo helps stabilizing the coalesced DNA 
stretch against reclosure, as the combined free energy of 
the two individual bubbles carries a factor CTq while the 
coalesced bubble has only a factor ctq- 

The full partition function is 

2f{X,Y) = 3fLiX)3fRiY). (3) 

It defines the free energy landscape J^{X,Y) = 
—f3~^\og[^{X,Y)], in which the random motion of 
the zipper forks takes place, as the gradient of ^ with 
respect to the coordinates X and Y defines the local 
driving forces experienced by the two zipper forks. 

Below the melting temperature of the barrier T^, the 
barrier will on average be driven towards closure, while 
above Tf, it will tend to denature completely. The effect of 
thermal fluctuations is to introduce a random walk-type 
dynamics of the position of the two zipper forks. Eventu- 
ally, full denaturation of the clamp may be reached even 
below the melting temperature T^. Once the two bub- 
bles coalesce, the loop initiation (cooperativity) factor 
(To is released, and the coalesced state becomes stabilized 
against closure. 

Dynamically we quantify the random motion of the 
two zipper forks due to thermal fluctuations as follows. 
To zip close an already opened base-pair, we assume that 



this process is mainly governed by diffusion-limited en- 
counter of the two separated bases, and subsequent bond 
formation. In contrast, to unzip a still closed base-pair, 
the free energy barrier embodied in the Boltzmann fac- 
tor u has to be overcome. For the left zipper fork we 
define < J (X, Y) which is the transfer coefficient for the 
process X X + 1, corresponding to clamp size de- 
crease, and t2 {X, Y) the transfer coefficient for the pro- 
cess X X — 1 (clamp size increase). For the right 
zipper fork we similarly introduce t~^{X,Y) for the pro- 
cess Y — ^ Y+1 (clamp size increase) and t~j^{X, Y) for the 
process Y ^ Y —1 (clamp size decrease). Due to the end 
clamping we require that X > —Nl and Y < N + Nr, 
which amounts to introducing reflecting boundary con- 
ditions [nl 

tl{X = -NL,Y) = 0, (4) 

and 

t+{X,Y^N + NR)^0. (5) 

Once the clamp has vanished, we assume that the clamp 
will not be able to reform for a long time, and we impose 
the absorbing conditions 

tl{X,X)=tUY,Y) = 0. (6) 

Dynamically, this is connected to the time it takes the 
long stretch of single-strand to re-establish a base-pair 
in the clamp region (diffusion limit). In terms of the 
free energy the suppression of clamp reformation is due 
to the release of the free energy AGbubbic corresponding 
to the cooperativity factor ctq on bubble coalescence (it 
would cost the additional factor (To to reintroduce two 
single-strand/double-strand boundaries) . 

Knowledge of the transfer coefficients together with the 
boundary conditions above completely determine the dy- 
namics, and we proceed by giving explicit expressions for 
the transfer coefficients in terms of the physical parame- 
ters of the problem. For the zipping rates we choose 

tl{X,Y) = ^IC{X + NL), (7) 

for the left fork, and identically for the right fork 

t+{X,Y)^^IC{NR + N-Y). (8) 

We defined above a bubble-sizc-depcndent rate coefficient 

IC{q) - fc(7-^ (9) 

with q being the number of broken bps in the bubble, 
where we have, as in previous studies, introduced the 
hook exponent ^, related to the fact that during the zip- 
ping process not only the base-pair at the zipper fork is 
moved, but also part of the single-strand is dragged or 
pushed along. This additional effect may be included us- 
ing similar arguments as in Refs. [2l|, |4l|, [s^ : To zip close 
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a base-pair, the two single-strands making up the bubble 
have to be pulled closer towards the zipper fork. The ad- 
justment of pulling propagates along the contour of the 
chain until the closest bend (inflexion) is reached, a dis- 
tance that scales as the gyration radius, i.e. ~ q'^ . Hav- 
ing in mind Rouse-type dynamics, this would slow down 
the unzipping rates by the factor q^'^. Hydrodynamic 
interactions may change the exponent and we here take 
the transfer coefHcients above proportional to g"*^, with 
/i to be determined by more detailed microscopic investi- 
gations. The rate constant k appearing in Eq. ^ is the 
rate constant for pure base pair unzipping without factors 
due to the coupling along the chain, i.e., the hook expo- 
nent. The factor 1/2 introduced above in Eqs. (O, ([S]) is 
merely for convenience to be consistent with the nomen- 
clature of previous approaches [2l|, |4l[ . Apart from the 
hook effect, we thus assign a factor k/2 for the zipping 
at each of the forks. 

As the DNA construct is embedded in a thermal bath, 
we require the zipping rates to fulfill the detailed balance 
conditions 

t+{X - 1, Y)^(X -!,¥) = tl{X)3r{X, Y), (10) 

and 

t],{X, Y + l)3f{X, Y + l) = t+{X, Y)fr{X, Y). (11) 

These conditions guarantee the relaxation to the ther- 
modynamic equilibrium. For the left zipping rates this is 
fulfilled for 

t+iX,Y) = ^ICiqL + lHX+l){{qL+l)/{qL+2)r. (12) 

Here qL ~ X + Ni is the length of the left bubble, and 
again {bubble length + is the correction for the en- 
tropy loss of a closed polymer loop, with c being the 
loop exponent. For bubble size increase we thus take the 
transfer coefficients to be proportional to the Arrhenius- 
factor u{X + 1), multiplied by a loop correction factor. 
Note that for q^ ^ oo the loop correction factor tends 
to 1, which we will exploit later. For the right fork we 
similarly find 



tRiX,Y) 



i/c(te + i)u(r 



2)r, 

(13) 

where qu = N + Nr — F is the length of the right bub- 
ble. We point out that Eqs. 0, dS]), and ^ are 
not unique in satisfying the detailed balance conditions, 
Eqs. (fTO|) and (fTTjl . However, different choices correspond 
to redefinitions of the time unit 1/fc which is the free pa- 
rameter in our master equation approach and needs to 
be fixed from fit to experiment. 

From the transition rates a master equation can be con- 
structed for the conditional probability P{X, Y; t\Xo, Yq) 
with Xq, Yq being the initial positions of the zipper forks. 
This master equation can be solved numerically, details 
of which are introduced in Sec. IVIl and physical quanti- 
ties such as the mean-first-passage time density can be 



calculated. However, based on four assumptions concern- 
ing the transition rates it is possible to derive a continu- 
ous Fokker-Planck equation approximating the full mas- 
ter equation description; this is done in Sec. IHII From 
the Fokker-Planck approach we then derive numerical re- 
sults for the coalescence time density, and both numerical 
and analytic expressions for the mean coalescence time 
and the probability density for the coalescence position 
in Sees. II VI and IVl These three following sections provide 
details and extensions of our previous short work (53 | . 



III. THE FOKKER-PLANCK APPROXIMATION 
TO THE MASTER EQUATION 

In this section we derive a Fokker-Planck approxima- 
tion to the master equation based on the following as- 
sumptions: 

(i) The temperature T is so high compared to Tg that 
the base pairs in the soft zones remain unzipped at 
all times, i.e. there are effectively reflecting bound- 
ary conditions at the interfaces between the soft 
zones and the barrier region; 

(a) the soft zones are sufficiently long such that the 
influence of the loop factors can be neglected, and 
similarly 

(Hi) the influence of the hook factors becomes suffi- 
ciently small. 

(iv) Finally, the number of bps in the barrier region is 
much bigger than one, i.e. iV 3> 1, which allows for 
taking the continuum limit (see below). 

Under the assumptions the full partition func- 

tion for the two bubbles and the partially denatured bar- 
rier region becomes [see Eq. ([3])] 



(14) 



where in this section X,Y G {0, 1,2,..., N} due to (z) 
with X < Y. Notice that X is the number of barrier 
base-pairs already broken from the left end of the barrier, 
and N — Y counts the broken barrier base-pairs from the 
right end. The free energy is given as 

^ = -2fcBr log <Jo - (Nl + Nr)AGs + {X + N-Y)AGb. 

(15) 

In the continuum limit (assumption (iv)) we introduce 
dimensionless coordinates x = ^ and y = j^, x,y & 
[0, 1]. The gradient of ^ with respect to the coordinates 
X and y defines the local force experienced by the two 
zipper forks, namely 



X 



dx 



= -NAGh, Fy 



d^ 
dy 



= NAGh = -F 



X, 



(16) 

and we immediately see, that the zipper forks X and 
Y are driven by opposite, constant forces as sketched in 
Fig. [3 
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Figure 2: Plot of the linear potentials experienced by the re- 
spective bubble interfaces in the case T < Tt (/ < 0) in terms 
of the dimensionless quantities x, y, f (see text for details). 

Using the simplifications (ii) and (Hi) stated above, the 
modified continuum rates (denoted by r to distinguish 
them from the notation introduced in the discussion of 



the discrete case) for closing a base-pair at the left fork 
at position x, or at the right fork at position y, become 
[see Eqs. dZ]), dH), lEl), and ^] 

rl{x,y)=r+ix,y) = k/2, (17) 

and 

rt{x,y) ^r~^{x,y) =Uik/2, (18) 

such that the zipping open of a base-pair requires crossing 
the barrier AGt,. The boundary conditions 

^1(0,2/) = r+{x,l) = r+{x,x) = r^(j/,y) = 0, (19) 

guarantee that base-pairs cannot close beyond the barrier 
region, and that the process ends when the two zipper 
forks coalesce. 

Define by P{x,y;T) the probability distribution that 
the left and right zipper forks are located at x and y, 
respectively, at some given time r. The time evolution of 
PCjg, ?^t) is then given in terms of the master equation 



-Q^P{x, y; r) = r+{x- 1/iV, y)P{x - 1/iV, y; r) + r" (a; + 1/iV, y)P{x + 1/iV, y; r) - [r" (x, y) +r+{x, y)] P{x,y- t) 
+r- (x, y + llN)P{x, y + 1/7V; r) + r+ (x, y - l/N)P{x, y-l/N-r)- [r^ (x, y) +r+{x, y)] P(x, y; r). (20) 

Following the standard derivation (sgI . [s?! we Taylor-expand the above master equation keeping the first two orders 
only. For instance, for the first term on the right hand side of Eq. (|20p . we obtain the Taylor expansion 

rtix^l/N,y)P{x-l/N,y;T) « r+(x, y)P(a;, y; r) - 1 A^+(x-, y)P(x, y; r) + ^-^r+ix,y)Pix,y;T). (21) 

This is the only consistent expansion of finite order according to the Pawula-Marcinkiewicz theorem (ssl . Issl . Is^ . 
Alternatively the full Kramers-Moyal expansion needs to be taken along. With analogous expansions for the other 
terms and after some rearrangement, we find the bivariate Fokker-Planck equation (ssj 

d f d d \ f \ 

— P(x, y; T\xo,yn) = F ( — - — I P{x, y; t\xo, yo) + D \ —-^ + ] P{x, y; r|a;o, yo), (22) 



r 



where, instead of the probability density P(x,y;T), we 
use explicitly the notation P{x, y; t\xo, yo) including the 
initial conditions Xq and yg. In Eq. (P^ . the force F and 
diffusion constant D are defined by 



and 



D 



kjuh - 1) 
2N 



^ k{ub + 1) 
47V2 



(23) 



(24) 



Eq. is completed by specifying the initial and bound- 
ary conditions. As initial condition, we choose the sharp 



(5- form 



P{x, y; 0|xo, yo) = S{x - xo)S{y - yo) 



(25) 



with xo < yo, and due to the initial condition Eq. (^5)1 the 
joint probability density P(x, y; t|xo, yo) is actually the 
Green's function of Eq. (P^ . The condition that the two 
bubbles in the soft zones are always open is guaranteed 
by the reflecting boundary conditions (here, we define 
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2/ = F/D) 



^ - 2/ ) P(a;,y;T|xo,?;o) 



= 0, 



x=0 



^ + 2/ ) P(a;,y;T|xo,?;o) 



= 0, 



(26) 



at the edges of the hne segment [0, 1]: Once a zipper fork 
reaches either edge, the only possible direction to move 
is to restart unzipping the barrier. Moreover, we specify 
the viciousness condition [73 



(27) 



according to which the two zipper forks cannot be at the 
same point: the two forks annihilate and the bubbles co- 
alesce. This last conditions ensures the continuous char- 
acter of the probability density P. For completeness, we 
actually need to specify a second set of boundary condi- 
tions. However, due to the viciousness condition (|27|) . we 
can choose this boundary condition ad libitum] a clever 
choice will turn out to be 



d_ 

dx 



P{x,y]T\xo,yo) 



— P{x,y;T\xo,yo) 
dy 



0, 



= 0. 



(28) 



Such a choice is possible because the zipper forks never 
reach these two points. 



IV. SOLUTION OF THE VICIOUS WALKER 
PROBLEM 

A. Transformation of the Fokker-Planck equation 

To obtain the solution of the Fokker-Planck equation 
(22), it is convenient to notice that after a redefinition 
of time unit t = Dt the problem depends on a single 
dimensionless parameter 



■' 2D Ub + 1 



(29) 



It is important that this parameter depends on the length 
of the barrier and the Boltzmann factor for opening 
the barrier bps but not on the kinetic constant k. Thus, 
apart from an overall prefactor fixing the time unit, the 
solution depends solely on the structural properties of 
the physical system under study. 

Let us summarize the rephrased problem in terms of / 
for completeness 



^ + + ^ _ 2/— + 2/— 
dt dx'^ dy'^ dx dy 



P{x,y;t\xo,yo) = 0, 
(30a) 



with boundary conditions 

d_ 
dx 
d_ 
dy 



2/ Pix,y;t\xo,yo) 



"^f j Pix,y;t\xo,yo) 
d 



^^P{x,y;t\xo,yo) 
d 

—P{x,y;t\xo,yo) 
dy 



0, 
0, 
0, 
0, 



y=0 



the viciousness condition 

P{x,x;t\xo,yo) = 0, 
and the ininitial condition 

P{x, y; 0\xo, yo) = 5{x - xo)S{y - yo) with xq < yo 



(30b) 
(30c) 
(30d) 
(30e) 

(30f) 
(30g) 



To proceed in the solution, let us first introduce the 
Fokker-Planck operators 



^ d'^ d 

^pp(a;) = ^-^ - 2/—, 



dx'^ 

d^ 



dx 
d 



Lpp(j.)-^ + 2/-^^, 



(31a) 
(31b) 



so that the Fokker-Planck equation (|30ap can be recast 
into the following form 



d_ 

dt' 



P{x,y;t\xo,yo) = h^p{x) +L,pp{y) P{x,y;t\xo,ya). 

(32) 

The operators Lpp(a;) and Lpp(j/) can now be used to 
transform our Fokker-Planck equation following the pro- 
cedure outlined in Ref. ^55!], Chapter 5.4. Let us now 
define the Hermitian operator 



.{x,xo) = e--^(^-^"'L+p(a;)e^(^-^«' 



d^ 

dx"^ 



-f. (33) 



Here the last equality sign can be shown by applying the 
operator to a test function. The operator L(x,a;o) corre- 
sponds to a Hamilton operator of a Schrodinger equation 
with imaginary time —iht, mass m = /2 and a constant 
potential The idea is that it is (in general) easier to 
solve the time-dependent Schrodinger equation than the 
Fokker-Planck equation, because the first-order deriva- 
tive has been eliminated. The relation between the solu- 
tions of the original and the transformed Fokker-Planck 
equations are found after a few lines of algebra, and we 
obtain the following result: li p^{x;t\xo) is a solution of 



d_ 

dt 



p^{x;t\xo) = l.{x,xo)p^{x;t\xo), 



then the density 



p^{x]t\xo) = exp[/(x- a;o)]p^(a;;i|xo) 



(34) 



(35) 
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solves the original equation 
d 



dt 



p (a;;i|xo) = Lpp(x)p {x;t\xo) 



(36) 



The boundary conditions in x, Eqs. (j30bp and (jSOdp . are 
transformed to 



^-f)p^ix;t\xo) 



dx 



+ /U^(x;t|xo) 



x=0 



x=l 



0, 



0, 



(37a) 
(37b) 



while the initial condition remains unchanged, p^{x;t = 
0\xo) = exp[-/(x - xo)]p^{x;t = 0\xo) = S(x- xq). In 
other words, with Eqs. ([33|l . (|37ap . and (j37b|l we have 
indeed transformed the original Fokker-Planck operator 
Lpp(a:) into Hermitian form. 
Noticing that 



L(y,2/o) 



f = ef^y~y'>'^i.^^{v)e-f^y-y"\ (38) 



the same procedure can be carried out for the y- 
coordinate, where the different sign in front of / in 
Eq. (|30ap causes that, if p^{y;t\yo) is a solution to 



^p^(2/; t\yo) = Hy, yo)p"{y; t\yo), 



then 



p"{y; t\yo) = cxp [-f{y - yo)]p^{y; t\yo) 
satisfies 

d 

— p^(2/;t|yo) = l^ppiy)p^{y;t\yo), 
and the boundary conditions are 

+ f)p''iy;t\yo) 
-Ap^'iy^tlyo) 



d_ 
dy 
d_ 
dy 



0, 
0, 



(39) 
(40) 
(41) 

(42a) 
(42b) 



y=o 



We now see why Eqs. pOd]) . pOe)) are clever choices for 
the additional boundary conditions on P(a;, y; i|xo, j/o)j 
together with the similarity transformations ([35]) and 
(|40p : reflecting the uneven symmetry of the problem with 
respect to the original coordinates, all equations defin- 
ing the functions p^{x] i|a;o) and p^{x; t\xo) are identical, 
and thus arc the functions themselves 



p^{x; t\xo) = P^{x; t\xo) = p{x; t\xQ). 



(43) 



B. Solution of the transformed Fokker-Planck 
equation 

After the individual similarity transformations for 
the left ([35]) and right (|40|) walker, respectively, are 
performed, we arrive at the following time-dependent 
Schrodinger equation 



d_ 

dt 



together with p^(y; t = Q\yo) = 5{y - yo). 



P{x, y; t\xo,yo) = [M^^ ^a) + My^ yo)] P{x, y; t\xa, yo), 

(44) 

with imaginary time; here, 

P{x, y; t\xo, yo) = e-/(--«)+/(^'-'/")p(a;, y; t\xo, yo), 

(45) 

which implies the same viciousness condition 
P{x,x;t\xo,yo) = as in the original formulation. 
Now, however, we effectively have two identical vicious 
walkers moving in a common (constant) potential for 
which the solution is well known [sol . [sil . It is given by 
the antisymmetric product of the single- walker solutions, 
namely 

P{x, y; t\xo,yo) = p{x; t\xo)p{y; t\yo)-p{y; t\xo)p{x; t\yo). 

(46) 

Note that this form of the solution is analogous to con- 
structing the solution of an absorbing boundary value 
problem for a single diffusor under a constant drift ac- 
cording to the method of images [g^I ■ 

The backward transformation of the solution (|46|) by 
inverting Eq. (|45p finally produces 



P(x, y; t\xo,yo) = e-^(-'-")^/(^-^") [p(x; t\xo)p{y; t\yo) - p(y; t\xo)pix; t\yo)] , (47) 

and by construction this is a solution of Eq. pOap . satisfying the boundary conditions Eqs. (j30bp - p0ep . as well as 
the viciousness condition P{x, x; t\xo, yo) ~ 0. It remains to check that this solution also satisfies the initial condition 
pOgP and, thus, fully solves the studied problem. For t = we obtain 

P{x, y; t = 0\xo, yo) - ef^^~-o)-fiy-yo) [^(^ _ ^^^s^y _ y^) „ s{y - xo)5{x - yo)] 
= Six - xo)6iy - yo) - e^f^y°-^°H{y - xo)5{x - yo) 

= 5{x-xo)5{y-yo), (48) 
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valid for xq < yg and x < y. The last equality follows 
because of the choice of the initial condition xq < y^ and 
the viciousness of the process, i.e. the walkers can never 
pass each other, which also implies x < y for any realiz- 
able configuration at any time t. Thus, the arguments of 
the (5-functions in the second term can never vanish si- 
multaneously and therefore this term is effectively equal 
to zero. Thus P7)) is the solution to our problem. 



C. Calculation of p{x;t\xo) and its spectral 
resolution 

To find the solution of the full 2-walker problem (|T7)) we 
need to calculate p{x; t\xo), which solves the Schrodinger 
equation 



9 ~, IN 

—p(x;t\xn) 



8^ 
dx'^ 



p{x;t\xo), 



(49) 



and satisfies the boundary conditions Eqs. (|37ap and 
P7b[) . as well as the initial condition p{x;t = 0\xq) = 
d{x — Xq). Unfortunately, there is no explicit solution 
of this equation in the time domain. It can be found, 
however, in the Laplace picture which is done in detail 
in Appendix |Al Here we only summarize the final result 



p{x;z\xo) = 



1 



2k{K^e^'' - I) 
+Ke 



Ke e 



2 2k ~k\x-xa\ 



with z-depcndent 



and 



k = k{z) = V^+p, 
k{z) + f 



K = k{z) 



k{z)-r 



(50) 



(51) 



(52) 



The corresponding behavior in the time domain is found 
by the inverse Laplace transform 



dz 

p{x;t\xo) = / —e''*p{x\z\xQ,) 

OO 

= Y^e^-'MxWni^o) (53) 



n=0 



with e > large enough such that all singularities of 
p{x; z\xo) lie in the half-plane Kez < e. We used the for- 
mal eigenmode expansion in the second line of the above 
equation which exists as a spectral resolution of the Hcr- 
mitian operator defined by Eq. (|49p and the pertinent 
boundary conditions (|37ap and (|37bp . The eigenvalues 
A„ are real numbers since the operator is Hermitian, but 
not necessarily non-positive like in the case of standard 
Fokkcr-Planck operators in one dimension. The reason 
is that p(x;t\xo) is an auxiliary mathematical quantity 
without any direct physical meaning and, thus, can in 
principle grow exponentially over time, i.e., some of the 
A„ might be positive. The physical quantity which is not 
allowed to grow indefinitely is the 2-walker Green's func- 
tion (j47p . This is useful to keep in mind when studying 
the spectral resolutions of fSOl) and (|47l) in more detail 
in the following. 

The spectrum of the operator d^/dx^ — from 
Eq. (P^ . together with the boundary conditions, can be 
found either directly from the defining equations or by 
finding the poles of the 1-walker Green's function (|50p . 
Indeed, the secular equation obtained by either method 
is equivalent to the denominator in Eq. (jSOp being zero, 
i.e., K^(A) ~ g-2fe(A)^ rpj^jg iga^jg transcendent equa- 
tions for real A in the respective ranges 



(A + 2f) sinh VA + /2 + 2/ v/A + p cosh y^xTp =0 for A > -/^ 
(A + 2f) sin v/|A| - p + 2/ VIAI - P cos ^lAI - P =0 for A < -/^ . 

I 



(54a) 
(54b) 



One looks for solutions of these equations in their 
ranges of validity. Thus, with the help of the standard 
graphical analysis, it turns out that the first equation 
(|54ap has at least one solution only for / < 0, this so- 
lution is positive, i.e. Ao > 0. When / < —2 a fur- 
ther solution with —p < Ai < solves Eq. (|54ap . 
There are no more options for Eq. (|54ap . On the other 
hand, the second equation (j54bp always has an infinite 
number of solutions. For / > these solutions are 



bounded by —f" 



1 



< A„ < 



-P 



{n — 0,1,2, . . . ). For > / > — 2 the first eigenvalue Aq 
satisfies Eq. (j54ap while the remaining eigenvalues (A„ for 
n = 1,2,...) stem from Eq. (|54bp and arc still bounded 
by ~p - n^TT^ < \n < -P - (n - 1)^tt^. Finally, for 
/ < — 2 the first two eigenvalues Aq and Ai are deter- 
mined by Eq. (|54ap and the rest stems from Eq. (|54bp 



being bounded by — / 
for n = 2, 3, . . . . 



<K<-P-in-iy 



There are two special values of the force / = 0, —2 
where the spectrum appears to change its analytic struc- 
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ture. First, the / = case corresponds to the problem 
of two vicious walkers freely diffusing in an impenetra- 
ble well, i.e. in a common potential. For this case, we 
do not need to use our trick as it is solved already by 
a previous study [sij : The single- walker spectrum reads 
An = (7i7r)^ for n = 0, 1, 2, . . . The second case, / = —2, 
does not appear to be in any way particular physically. 
Actually, neither of the two cases have any exceptional 
physical properties — all the physical quantities change 
smoothly across these two points when changing /. The 
apparent singularities occur only in the auxiliary quanti- 
ties. Technically, these two cases are the only ones where 
the trivial solution A = — of Eqs. ([54)) corresponds 
to a non-trivial, i.e. non-zero, solution for the eigenfunc- 
tion. Mathematically, this is reflected by the formal fail- 
ure when k{z) = of the method used in Appendix |^ 
leading to Eq. ([50]) . In such cases, the two fundamental 
solutions e^'"^ of Eq. (jA3[) become identical and equal 
to a constant. The second independent solution is then 
linear in x according to the elementary theory of linear 
differential equations. Thus, the two special cases are 
not directly covered by the general solution presented in 
App. [X] and all formulas stemming from it. We do not 
give explicit solutions for those two singular cases since 
all quantities of interest can be obtained from the general 
formulas by taking the appropriate limit / — > 0, or —2. 

Now, if we insert the expansion ([55]) into Eq. ((Tf|) 
we see that the spectrum of that equation is given by 
all pairwise sums of different 1-walker eigenvalues, i.e. 
Afc = Xi+Xj {i 7^ j) [the sums of identical 1-walker eigen- 
values have zero weight due to the anti-symmetrization 
in (|47)) ]. Therefore, the first two 1-walker solutions 
Aq and Ai determine the long-time asymptotics of 
the 2-walker problem. For large negative / ^ —1, 
corresponding to a large barrier case, the asymptotic 
behavior is dominating the whole solution and, thus, 
the combination Aq = Aq -I- Ai effectively determines 
the interesting quantities such as the mean coalescence 
time or the probability distribution of the coalescence 
position. This allows us to get full analytic results in 
the limit / ^ — 1 , as explicitly derived in Section IV Bl 
In this case both eigenvalues Ao,i are solutions of 
Eq. (|54ap . and they are exponentially small in |/|. They 



can be found to leading order from the second order 
expansion in A of Eq. (|54ap . yielding Ao,i — ±4/^6^1-^1. 
To obtain Aq we need to increase the accuracy by 
expanding the equation up to the third order in A and 
we find Ao = Ao + Ai ~ ~lQp{\f \ - l)e-2|/l which is 
negative, as it should be since P(a:, y; tjxo, yo) cannot 
grow exponentially, P being a physical quantity. So, 
despite the existence of a positive 1-walker eigenvalue 
Ao, the physically relevant combination Aq + Ai < 
ensures meaningful results. 

V. BUBBLE COALESCENCE TIME AND 
POSITION FROM THE SOLUTION OF THE 
FOKKER-PLANCK EQUATION 

A. General treatment 



From the solution of the Fokker-Planck equation, we 
now calculate in general the quantities of interest, namely 
the characteristic coalescence time of the two random 
walking zipper forks, and the probability distribution of 
the coalescence position. In the time domain this is done 
by considering the conservation of probability in the form 



n(t|a;o,?/o) + / dy I dxP{x,y;t\xo,yo) = (55) 



which is actually a defining equation for the probability 
n(t|a;o, yo) that the walkers have met before time t. Con- 
sequently, the second term represents the probability of 
having two separate bubbles at time t (survival probabil- 
ity). Note the range of integration restricting x £ [0,y)- 
The probability density associated with Il(i:]xo,yo) is 



7r(i|a;o,yo) = -jj'n{t\xa,ya) 
dt 



1 ry g 
dy / dx—P{x,y;t\xo,yo), 



(56) 



which, after invoking Eq. (|30ap and rearranging the order 
of the integrals, yields 



7r(i|xo,?;o) = - 



1 ry 
dy / dx 

"'0 

1 



dx'^ dx 



P{x,y;t\xo,yo) - / dx dy 
Jo Jx 



^ + 2/ A 
dy'^ dy 



dy 







d 

— - 2/ \ P{x,y;t\xo,yo) 



"'O 







I ^ - 2/ J- P{x, y\ t\xQ, yo) 



dx 1^-^ + 2/} P{x,y;t\xo,yo) 
dx I S- + 2f'> P{x, y; t\xQ,yo) 



P{x,y;t\xo,yo) 
1 







dy 



(57) 



y=x 



,yo) 



y=x 
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where the boundary conditions imposed on 
P{x,y;t\xo,yo) and the viciousness condition have 
been used in subsequent manipulations. The quantity 
7r(<|a::o, yo) can be interpreted as the probabihty current 
flowing into the absorbing boundary, here given by the 
hne X ^ y (compare with an analogous discussion for a 
1-dimensional case in Ref. [6l|V 

From the last relation, we define the quantity 
g{x;t\xo,yo), the coalescence-position-resolved probabil- 
ity density for the coalescence time, 



Tr{t\xo,yo)^ / dx g{x;t\xo,yo), 



(58) 



such that 



g{x;t\xo,yo) 



d_ 

dy 



d_ 

dx 



P{x,y;t\xo,yo) 



(59) 



We are mainly interested in either the mean coales- 
cence time T{xo,yo) = dttn{t\xo,yo) or the proba- 
bility density of the coalescence position p{x\xo,yo) = 
dt q{x; t\xo,yo)- These are quantities integrated over 
time and, thus, they may be determined from the solu- 
tion in the Laplace domain without explicit knowledge of 
P{x, y; t\xo, yo) in the time domain. 

We now rewrite Eq. (|47p in terms of the inverse Laplace 
transforms as follows 



pt+ioo 1 pe+ioo j 

Pix,y;t\xo,yo) = 77^^^'' -^e^^*e/(--^--''+«'') [p(x; zi|xo)p(y; z^lyo) - p(y; zi|xo)p(x; |yo)], (60) 

where p{x]z\xq) is given in ()50|) . In case that p{x]z\xq) has no poles in the half-plane Rez > (corresponding to 
/ > 0, e = 0), we can directly use the substitution zi = z/2 + icu, Z2 = z/2 — ilu; z = zi + Z2, and uj = (zi — Z2)/(2i), 
to obtain a single inverse Laplace transform for P{x, y; t\xQ, yo), namely 

/ioO 7 pOO 7 

_^^fi.-y-.o+yo) [p(:c;z/2 + tu;\xo)p{y;z/2-iLu\yo) 

-p {y; z /2 + iuj\xo) p {x; z /2 - iuj\yo)] . (61) 
Thus, we can identify the Laplace transform P{x, y; z\xq, yo) of P{x, y; t\xQ, yo) for / > as 

JH^fi^-y-'-^o+yo) [^(2.. + iuj\xo)p{y; z/2 - iuj\yo) ^ p{y; z/2 + iuj\xo)p{x; z/2 - iuj\yo)]. (62) 



This directly yields, by means of relations ((58l) and ((59)) . 
the Laplace transforms tt (zjxo, j/o) and g{x; z\xo,yo), 
from which we in turn deduce the (time-averaged) dis- 
tribution of coalescence positions, 

p{x\xo,yo) ^ g{x;z ^0+\xa,ya), (63) 
and the characteristic (mean) coalescence time 



T{xo,yo) = - ^7r(z|xo,yo) 



(64) 



z=0+ 



For / < the situation is more complicated since 
p{x; z\xq) in this case has a pole at Aq > which prohibits 



us from just repeating the above reasoning. The Laplace 
transform (|60p only holds for Rezi^2 > Aq > implying 
Re z > 2Ao > and the analytic continuation down to 
z = is not obvious. However, since we know the posi- 
tion of the only pole located in the Re z > half-plane, 
we can treat this singularity separately and thus gen- 
eralize the previous results. Using the Cauchy theorem 
for complex integrals we move the integration line from 
€ + iuj with e > Ao > down to the imaginary axis and 
add the contribution from the (single) singularity at Aq in 
between these two lines. This method eventually leads to 
an expression for the Laplace transform P{x,y; z\xo,yo) 
in the whole half-plane Rez > for / < 0, 
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C /"OO T 

P{x, y- z\xo, yo) = e/(--f--"+««) | j z/2 + *c^|xo)p(y; z/2 - ic^|yo) - z/2 + *u;|xo)p(x; z/2 - iio\yo)] 

+ Rcs2i=AoP(a;; 2i|a;o)p(y; 2: - Ao|yo) - Rcs2i=AoP(y; zi|xo)p(a;;z - A0I2/0) 
+ p(a;;z - Ao|a;o)ReSz2=AoP(y; 22|yo) -p{y;z~ Xn\xn)ReSz2=XoP{x; Z2|2/o)|, 

(65) 

which is the sought-for generahzation of ([S^ . We observe that the double- residue term drops out due to the anti- 
symmetrization procedure. 



r 



B. Results 



Eqs. (|62l) and ^ for / > and / < 0, respectively, 
together with the single- walker Green's function ([50]) and 
the identities for p{x\xo,yo) (j63p and T{xo,yo) (j64p were 
implemented in MATHEMATICA and evaluated. The results 
are shown in Figs. [3] and H] depicting the coalescence po- 
sition probability density for several values of / and the 
mean coalescence time function of /, respectively. 
Both quantities are shown for two different initial con- 
ditions xq = 0, yo ~ 1 (the two walkers start out right 
at the boundaries) and xq = 0.5, yo = 0.9 (a generic ini- 
tial condition). Further results for the coalescence time 
probability density 7r(i|a;o, yo) (|57p are presented in later 
sections. 

The results for the coalescence position in Fig. [3] show 
for both initial conditions a clear crossover from a peaked 
form of the probability density for large positive force 
(the case of an almost "free fall" into the potential 
well where the boundary conditions have negligible in- 
fluence on the dynamics, studied in detail in Ref. [3l| 
and discussed below) to a very flat probability density in 
the case of large negative force corresponding to a high 
barrier. The flatness in the latter case can be under- 
stood from a simple Arrhenius-like model, in which the 
probability of the walker to be at a place x is propor- 
tional to a Boltzmann weight, i.e., exp[—f34>{x)], where 
(j){x) = — F{x')dx' is the free energy corresponding to 
the force F{x) = ±f. Since we are now dealing with two 
walkers, the probability of both of them being simulta- 
neously at the coalescence position is given by the prod- 
uct of the Boltzmann weights, cxp{~(3[<j>L{x) + <j>]i{x)]), 
which is a position-independent constant due to the can- 
celation of the position-dependence of the two opposite 
linear potentials. This simple picture breaks down close 
to the boundaries but otherwise is sufficient to grasp the 
observed behavior. 

The characteristic (mean) coalescence times in Fig. 2] 
cross over from the "free fall" behavior for / ^ 1, propor- 
tional to the inverse of the force r ~ 1// (using the natu- 
ral boundary condition in the calculation of Ref. [3l| gives 
''' — {yo — a;o)/(4/), cf. Eq. (9) therein) to the thermal 
Arrhenius/Kramers-like barrier crossing proportional to 
the exponential of the barrier height r ~ exp(2|/|) for 



/ — 1. Thus, all the results are plausible and can be 
qualitatively rationalized based on simple physical argu- 
ments. 

In the rest of this section we focus on a more detailed 
study of the two limiting cases, the almost "free fall" 
/ ^ 1 and the large barrier / ^ — 1. For these limit- 
ing cases we obtain analytical results from relatively sim- 
ple assumptions which compare quantitatively well with 
the full solution. We start with the "free fall" (ff) case 
where we assume that the large drift towards coalescence 
dominates the dynamics so that the reflecting bound- 
ary conditions can be safely neglected since the typical 
realizations/trajectories of the stochastic process never 
reach them. The validity of this assumption depends on 
the initial conditions and we expect it to be good enough 
for the walkers starting from initial positions xq-, yo sat- 
isfying min(a:o, 1 — yo) ^ far enough from the 
boundaries. This, indeed, turns out to be the case, see 
below and compare the corresponding results in Figs. [3^ 
and[3jD. 

If the reflecting boundary conditions are neglected the 
2-walker Fokker-Planck equation (|30a[) can be solved by 
separation of variables. In particular, the reformulation 
of Eq. (|30ap together with conditions (jSOfp and pOgP in 
terms of the center-of-mass {x + y)/2 and relative y — x 
variables leads to a separable problem which can be easily 
solved since the center-of-mass coordinate (cms) just per- 
forms a free diffusion while the relative coordinate (rel) 
satisfies equations analogous to those in Ref. [3l| (Eq. 
(7) with c=0 therein). Following a derivation similar to 
that in Sec. IV Al we arrive at (for details the reader is 
referred to an upcoming publication [6^ ) 

Qfi{x]t\xo,yo) = -Terns ( ^ ■,t\ — ^ — jTrroUqyo " a;o), 

(66) 



with 



-Perns (w; t\ua) 



: CXp 



[u - Up)' 

2i 



(67) 



being the free diffusion propagator of the center-of-mass 
coordinate (see Refs. [63,|6J|) and 



X yo - 2^0 

7rrcl(i|yo ~ Xq) ^ CXp 



(yo - xq- Aftf 
8t 



(68) 
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Figure 3: Probability density for the coalescence position 
p{x\xo, yo) as a function of the position x for several values of 
the dimensionless force /. The initial positions of the walkers 
xo and yo were xo = 0, yo ~ 1 (a) and xo — 0.5, yo = 0.9 
(b). The full lines for cases / = 10 and / = — 10 correspond 
to analytical results in Eqs. (|69|) and (|74|) . respectively, for 
given initial conditions. 



the first passage time probability density for the rela- 
tive coordinate to reach the origin (compare with Eq. (8) 
in Ref. [3l|). Using p{x\xo,yo) = dt g{x]t\xo,yo) 
and the identity dtexp{-a'^ / (2t) - 2b'^t)/t^ = 
46Ki(2a5)/a for a, 6 > (K„(x) is the modified Bcsscl 
function of the second kind of order n) we finally ob- 
tain for the coalescence position probability density in 



10= 









— x„=0, y„=1 








...x„=0.5,y„=0.9 








1/(4f) 








---0.4/(4f) 








^exp(2|f|)/[16f^(|f|-1)] 











0-1 . . . 

-10 -5 5 10 



f 

Figure 4: The mean coalescence time T{xo,yo) as a function 
of the dimensionless force / for two different initial conditions 
Xo = Q, yo = 1 (full line) and xo ~ 0.5, yo ~ 0.9 (dashed line). 
The asymptotic analytical results for large positive, i.e. "free 
fall" case with r — {yo — x-o)/(4/), and negative, i.e. large 
barrier case of Eq. f75)l , forces are also shown for comparison. 

the "free fall" limit 

/ I s ./(2/o-a::o)exp[/(?;o-a:^o)]Ki(2/r) 

n r 

with 

r = y/{x - {xo + yo)/2)2 + {yo - xof/^. 

(69) 

This function is plotted in Fig. [3|i,b for the two differ- 
ent initial conditions for force / = 10. Wc can see a 
rather good agreement between the asymptotic formula 
([55)1 and the full result in Fig. The situation is much 
worse in Fig. [3^ although even there the correspondence 
is qualitatively quite acceptable. As mentioned above, 
the reason for the success or failure of the approximation 
is determined by the initial conditions. Indeed, Fig. [3)3 
corresponds to /min(a;o, 1 — yo) = /(I — yo) = 1 where 
the approximation is expected to become valid while in 
Fig. [3^ the walkers start out right at the boundaries and 
only an extremely high value of the force could prohibit 
the walkers from occasionally bumping into the bound- 
aries, especially at the very beginning. Thus, in the case 
with xq = 0, ?/o = 1 the above approximation is expected 
to become quantitatively accurate only for very high val- 
ues of / — numerical estimates reveal that an agreement 
comparable with that of Fig. is not achieved until 
about / > 40. This supports a heuristic guess that the 
accuracy of the asymptotic analytic expression crosses 
over from exp[— /min(a;o, 1 — ?/o)] to 1// when the mini- 
mum equals zero. 

Now, we turn to the opposite limit of large barrier. 
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i.e. the / <C — 1 case. In particular, we want to derive 
asymptotic expressions for the characteristic coalescence 
time (which is identified from the results of the full the- 
ory as T = 1/1 AqI, cf. the end of Sec. lIVCl and below) and 
the coalescence position probability density. Clearly, this 
limit is dominated by the lowest eigenvalue and eigen- 
function of the 2-walkcr problem which means by the two 
lowest eigenvalues and eigcnfunctions of the auxiliary 1- 
walker problem. Thus we can write, under assumptions 
of large barrier (lb) and generic initial conditions (to be 
specified in more detail below), for Pih{x, y; t\xo, yo), 

Pib{x,y;t\xo,yo) ^ e^(^-v)erf'^'">-y°^e^^'>+^'^' 
X [■iPo{x)ipi{y) - ipo{y)ipi{x)] 
y-[tpo{xo)ipi{yQ) - tpo{yo)il!i{xo)],{7Q) 

where Ao,i — ±4f^e~^^^ are the lowest eigenvalues satis- 
fying Eq. (|54a|) and ipoAix) are the corresponding eigcn- 
functions given by (these are exact expressions for all 
/<-2) 



Tplix) 



2Ao 



Ao + 2|/| 

I 2|Ai| 
Ai+2|/| 



cosh[v7^TA;;(x-i/2)], 



(71) 



sinh[v//2 + Ai(a; - 1/2)] . 



Using these expressions we can study the dependence on 
the initial conditions and clarify the regime in which the 
assumption about the dominance of the lowest eigenmode 
is valid. Utilizing that in the limit / <C — 1 the eigenval- 
ues satisfy |Ai| ~ Aq <C |/| we obtain 

~ I/I 



^(1 

2i/r 



expl 



sinh[|/|(?;o - a;o)] 
-2|/|(yo-2:o)]) . (72) 



Since by assumption yo > xq we see that the evolution de- 
pends only exponentially weakly on the initial conditions 
so that for yo — xq ^ 1/|/| the evolution is essentially 
independent of the initial conditions as expected in the 
high barrier limit. Indeed, the above condition just says 
that the walkers should start out well separated so that 
the barrier between them is still large (in dimensionless 
units) . In such a case the dynamics is independent of the 
detailed initial condition or, more precisely, it depends 
on it only exponentially weakly which can safely be ne- 
glected. This is the regime in which the lowest eigenmode 
theory of Eq. ([70)1 is sufficient as we will demonstrate be- 
low. 

If we calculate g{x; t\xo, yo) from Eq. (j59|) in the limit 
/ < -1, yo - xo > l/l/l we obtain 

gii,{x;t\xo,yo) ~ ^e^^'+^'^'iiJoWiix) - M^Xi^)] 



X2 



I/I 



/I 

o + Ai)t| 



1 



2/2 



cosh[2|/|(x-l/2)]}. 



Now integrating over time and taking into account that 
Ao ~ 4pe-\f\ and Aq + Ai =^o =^ -16/2(|/| -l)e-2|/l < 
we get for p{x\xo,yo) = J^^ dt g{x; t\xo , yo) 



Pib(a;|xo,yo) 



I/I 



l/l-l 
1 



1 



I/I 



1 -2e"l-''lcosh[2|/|(2;- 1/2)] 



l_e-2|/l-_e-2|/l(i--)). (74) 



(73) 



This clarifies that pih{x\xo,yo) is properly normafized to 
one within exponential precision, dx pi\:,{x\xo,yo) = 
I + O (e^^l-'^l) ~ 1 , which finally proves the self- 
consistency of the lowest eigenmode approximation. The 
curves for the case / = —10 in Fig. [3^,b calculated by 
the full theory are practically indistinguishable from that 
given by the approximate expression (|74p (which is an ex- 
plicit illustration of the initial-condition independence). 
In a straightforward manner it also follows that the mean 
coalescence time is given by the inverse lowest eigenvalue 
1/1 AqI since for gii,(x;t\xQ,yo) in the separable form of 
Eq. (|73p and due to the above normalization condition 
one immediately gets 

Tib(a;o,yo) - / dtt dx gi^,{x;t\xo,yo) 
Jo Jo 

= / dtte"l^«l*|Ao| / dxpn,{x\xo,yo) 
Jo Jo 
1 e^\f\ 

^ JjQ^ 16/2(1/1-1)' ^^^^ 

independent of the initial conditions. All approximate 
equalities hold up to exponentially small corrections of 
order e"^'^', which are negligible for / <C —1. 



C. Summary 

We have in Secs. lIIIillVi andlVlset up an approximate 
Fokkcr-Planck equation scheme for the full problem of 
two interfaces moving in a block DNA-stretch with a 
barrier region separating two soft zones. While the full 
problem can be viewed as two discrete random walkers 
in different potentials with an imposed vicious boundary 
condition, the approximate Fokker-Planck equation de- 
scribes two continuous random walkers in opposite linear 
potentials keeping the imposed viciousness condition. 
The four assumptions leading to the Fokker-Planck 
equation were introduced in Sec. IIIII and their validity 
will be discussed favorably in Sec. IVIl] when we compare 
the Fokkcr-Planck results with the direct evaluation 
of the full discrete problem using the master equation 
approach presented in the next section. 

The main outcome of the Fokker-Planck approach are 
the general results for the coalescence time density 7r(t) 
(examples will be shown in Sec. IVIlj) . and the numerical 
and analytic expressions for the mean meeting time r 
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(shown in Fig. [4|) as well as the probability density p 
for the meeting position (shown in Fig. [3]). All results 
are expressed through the dimensionless force /, which 
in terms of the parameters from an experimental setup 
reads 



N{ut, - 1) 



Ub 



1 



(76) 



making comparison with values obtained from experi- 
ments straightforward. 

Finally, a remark on why one should consider the 
continuous approach over the complete, discrete mas- 
ter equation approach is in order. Namely, for DNA- 
stretches of length N the discrete master equation ap- 
proach involves diagonalization of matrices of the order 
N"^ X N"^ , setting computational limitations on N and 
a new diagonalization is needed for each parameter set. 
The Fokker-Planck approach may therefore provide addi- 
tional insight for very long DNA stretches as we showed 
here by discussing the physical quantities of meeting po- 
sition and meeting time. 



VI. COMPLETE DISCRETE APPROACH: THE 
MASTER EQUATION 

In this section, we develop a master equation frame- 
work for the bubble coalescence. In contrast to the pre- 
vious treatment, we explicitly allow the soft zones to zip 
close from the two ends of the barrier region. It will turn 
out that in some cases this only has a minor effect. A 
detailed comparison with the continuum Fokker-Planck 
equation approximation is shown in the next section. 

We consider the same segment of double-stranded 
DNA with M = Nl + N + Nr internal base-pairs, 
clamped open at both ends according to Fig. [T] How- 
ever, in contrast to the approximations imposed in the 
Fokker-Planck approximation, we now allow for explicit 
closure of the soft regions, necessitating the consideration 
of a sequence-dependence of the local DNA stability, in 
contrast to the previous discussion, where we assumed 
that the two bubble domains are always remaining open. 

Note that in this section our notation differs from the 
scheme introduced above, first, in order to keep the no- 
tation of this section consistent with previous references 
on the same method [2l|, [H, [H, [H, [sl] , and second, to 
be able to incorporate zipping/unzipping of base pairs 
in the soft zones, as well. We denote hy xl = X + Nl 
{xr = Y -\- Nl) the position of the rightmost (left-most) 
open basepair in the left (right) open region, see Fig. [U 
where xl,r G [0, M]. The positions xl and xr of the 
two zipper forks are stochastic variables and the aim is 
to understand how these variables evolve in time without 
taking the continuum limit and using the approximations 
introduced in Sec. IIIII We note that an equivalent set of 
variables are xl and the clamp size m, that are related 
through 





Axl 




Am 


Eq 


f^ixL, m) 


1 





-1 


(|8()|) 







-1 


-1 




tl{xL,m) 


-1 





1 


(|82|) 


t%{xL,m) 





1 


1 


(|83|) 



Table I: Properties of the transfer coefficients. The quantity 
Axl (Axr) denotes the change in fork position xl (xr) under 
the action of the transfer coefficients. Similarly Am denotes 
the change in clamp size m. 



In the master equation formulation below we will use x l 
and m as the dynamic variables, and for completeness 
we state the transition rates, Eqs. ([3]), ([5]), ([7]), (O, and 
(fT2|) . (fT3|) expressed in the new variables: The reflecting 
boundary conditions are 



tj^ {xl = 0, m) = t^{xL,rn 



M-xl) 



0. 



(78) 



Once the clamp is completely unzipped, i.e. the state 
TO = is reached, we assume that the clamp will not 
be able to reform for a long time, and we impose the 
absorbing conditions 

tl{xL,m^O) ^t+{xL,m^O) ^0. (79) 

The transition rates at the interior of the DNA stretch is 
1. 

2' 
1 

2' 
1 

2' 
I' 

xs{M — m — Xl), 



t^{xL,m) 



tl{xL,rn) 



tJxL^m) 



-IC{M - m - xl), 

^-JC{xl + l)u{xL + l)s(m), 

-JC{M — m — + \)u{xb 



1) 



(80) 
(81) 
(82) 

(83) 



where Xr = Xl + m -\- \, s{q) ~ {{q + l)/{q + 2)Y and 
JC{q) = kq~^. The properties of the four transfer coeffi- 
cients above are summarized in Table HI 

Denote by P{xl, m, t\x'j^, m') the conditional probabil- 
ity to find the system in state {xl,iti) at time t, given 
the initial condition (x'^,to') at time t = 0. With the 
short-hand notation P{xl, 'ti, t) = P{xl, m, t\x'j^, m') the 
dynamics is described by the master equation 



d_ 

di 



P{xL,m,t) = 

i+(.TL - 1,?71+ 1)P(.TL - 1,TO+ l,i) 
+tl{xL + 1,TO- 1)P{XL + 1,TO- l,i) 

- [t J [xl ,m) +tl (xl , m)] P{xL,m,t) 
+t+{xL,m - l)P{xL,m - l,t) 
+t]^{xL,m+ l)P{xL,m+ l,t) 

- [t'^{xL,m) +t]^{xL,m)] P{xL,m,t). 



(84) 



xr - XL- 



(77) 



This equation states that the probability for the clamp 
size can change in 8 different ways: the terms with a 
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plus-sign correspond to jumps to the state {xL,m}^ and 
the terms with a minus-sign correspond to jumps from 
the state {xL,m}. 

A standard approach to the master equation ([84|) is 
the spectral decomposition [H, 



A. Coalescence time density 

The (survival) probability that the absorbing bound- 
ary at m = has not yet been reached up to time t is 



P{xL,m,t) ^^Cp{x'i^,m')Qp{xL,m)exp{-T]pt), (85) {^'l , m' , t) = Pi^L,m,t\x'L,m'). (91) 



m— 1 xt—0 



in terms of the eigenvalues rjp and eigenvectors 
Qp{xL,rn). The expansion coefficients Cp{x'i^,m') are ob- 
tained from the initial condition. As in the previous sec- 
tion, we will assume the system initially to be in the state 
where all bps in the soft zone are broken and all bps in 
the barrier region are closed. 

The eigenvalue equation corresponding to Eq. 
becomes 

tJ(xL - l,?7i+ l)Qp{xL - l,m,t) 
+tl{xL + l,m- l)Qp{xL + l,m- l,t) 

- [iJ(a;L,m) -f t£ (ccl, 771)] Qp{xL,m,t) 
+t^{xL,m - l)Qp{xL, m- l,t) 
+t~^{xL,m+ \)Qp{xL,m + l,t) 

- [t%,{xL,m) + i^(xL, 77i)] Qp{xL,m, t) 

= -VpQp{xL,m,t). (86) 

The eigenvectors satisfy the orthogonality relation [s^ 

^'^^ Qp{xL,m)Qp'{xL,m) _ 

pcq^^. _A -<Jpp', (St) 



m—l xl—0 



Pr{xL,m) 



where 



with 



P;i(a;L,m) = 3f{xL,m)/2f, 



The probability that the absorbing boundary is 
reached within the time interval [t, t+dt] (namely, the co- 
alescence time density corresponding to the first passage 
problem) is 

p{x'L,m\t)dt ^ .y{x'i^,m',t)-y{x'L,m',t + dt) 

= - (^^^(x'i,m',t)^dt. (92) 

This expression is positive, as .y is decreasing with time. 
To express the coalescence time density in terms of the 
eigenvalues and eigenfunctions, we introduce the eigen- 
mode expansion into equation ([HH), yielding 



with coefficients 



Cp{x'L,m') 



VpCp{x'L,m')exp{-r]pt), 



(93) 



]p{x'^,m') 



M M-m 



TtE E Qpi^L,m). (94) 



(89) 



We have above made use of the orthonormality relation 
(|87p in order to express Cp{x']^,m') in terms of the ini- 
tial probability density P{xl, tti, Ojx'^^, m'), and used the 
fact that this general initial condition takes the explicit 
form P{xL,m,Q\x']^,m') = 5x^^x'^5m,m' ■ Eq. is the 
discrete counterpart of the continuous result derived in 
Ref. [131 , and expresses the coalescence time density (for 
any given initial condition, specified by m' and x'j^) in 
terms of the eigenvalues and eigenvectors of Eq. 



being the partition coefhcient in the variable xl and ?7i 
(neglecting a common cTg-factor, where uq is the bubble 
initiation parameter), and where 



Af-l-l 



u{xr) 

xl=0 XR=m+XL + l 

xil + XLy'^iM -m-XL + iy, (90) 



see Eqs. ©-(IS]). From the eigenvalues and eigenvectors 
of Eq. (|86p any quantity of interest may be constructed. 
In the following subsection we calculate the coalescence 
time density. How to set up the master equation for 
numerical purposes is presented in detail in App. [B] Al- 
ternatively, the master equation ([M)) can be solved by 
direct stochastic simulations such as, e.g. the Gillespie 
algorithm introduced in App. [C] and used in Sec. IVIII 



VII. COMPARISON BETWEEN THE FULL 
MASTER EQUATION AND THE 
FOKKER-PLANCK APPROXIMATION 

In this section, we investigate the validity of the as- 
sumptions presented in Sec. IIIII leading to the Fokker- 
Planck continuum approximation. This is done by com- 
paring the results for the coalescence time densities, 7r(t), 
with the results obtained from the full discrete master 
equation approach. In the next section, Sec. IVIIIl we 
discuss the relevance for biological experiments. 

In all examples below the two walkers move between 
[0, N] in the Fokker-Planck description and between 
[—Nl, N + Nf{] in the master equation setup. That is, in 
the master equation approach we explicitly allow zipping 
of base pairs in the two soft zones. As initial conditions 
we use Xq = and Yq = N throughout this section. 
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Figure 5: Coalescence time probability density for barrier 
width A'^ = 20 and temperature below and above Tt, with 
Uf, = 0.98 and ui, = 1.1, respectively. To exclude other effects, 
the lengths of the soft zones are zero, A''^ = Nr = 0, i.e., these 
are assumed to be always open, and c = fj, = 0. 

A. The continuum approximation 

The continuum assumption (iv) implies that the in- 
herently discrete nature of the DNA structure - both in 
terms of stacking and hydrogen bonds - can be approx- 
imated by the diffusive behavior of two continuous vari- 
ables. To get the Fokker-Planck description one has to 
consider the limit a — > with a being the length between 
effective bonds in the base pairs. In practice this limit is 
obtained by 

bond distance ^ 
total segment length ' 

i.e. by increasing the width of the barrier region. The 
question is addressed in a setup with open soft zones, i.e. 
Us ^ 1, and varying barrier lengths, A^. To have per- 
fectly reflecting boundary conditions we set A^^ ~ Nji = 
in the master equation setup. 

Fig. [S] shows that it requires a relatively small number 
of base pairs (~ 20) before the continuum approximation 
is reasonable, independently of the temperature. Barrier 
regions of this length arc in principle accessible experi- 
mentally so the continuum approximation appears to be 
well-justified. 

B. Open soft zones 

The first assumption, (i), states that the soft zones are 
always open, i.e. that the random walkers are reflected 
at the interfaces between the barrier region and the soft 
zones. 

To eliminate other effects than the effect of the intro- 
duction of the reflecting boundary conditions at the ends 



of the barrier region, we consider a DNA stretch of length 
25, so that the continuum approximation is justified, and 
set c = /.t = in order to exclude effects originating from 
the entropy factor and the hook exponent. We compare 
this to the results from the full master equation including 
the soft zones. 

Fig.[6]shows the coalescence time density 7r(i) for vary- 
ing lengths of the soft zones for temperatures above and 
below the melting temperature of the barrier region. Ap- 
parently, the length dependence is rather weak as long as 
the soft zones serve as hard enough boundaries, i.e. for 
large enough Ug > 5. For smaller Us the soft-zone- length 
dependence is relevant as the two bubble corners ven- 
ture more frequently into the soft zones. This, however, 
implies the breaking of the assumptions for the applica- 
bility of our Fokker-Planck description as revealed in the 
figure. Thus, the length of the soft zones itself is not im- 
portant if Us is large enough and c = fi = 0. Systematic 
assessment of these conditions is given below. 

Fig.[7]shows TT{t) as a function of Ug for A''^ = = 20, 
i.e. the soft zones are so long the two forks essentially 
never reach the outer clamps. That this is indeed the 
case can be qualitatively investigated using the Gillespie- 
scheme presented in App. [Cl giving access to real-time 
trajectories of the two random walkers in a potential 
landscape including both the barrier region and the soft 
zones. This is shown in Fig. [5] which confirms that ex- 
cursions into the soft zones are progressively suppressed 
with increasing Us- In Fig. [7] we have only included 
the case Ub < 1, i.e., when the barrier region indeed 
acts as a barrier, and consequently the effect of the soft 
zones is more pronounced. The discrepancies between 
the Fokker-Planck and the master equation approach be- 
come distinct for ~ 1 whereas for > 5 the agree- 
ment between the approaches becomes reasonable. Dif- 
ference between Us and Ub of this magnitude can indeed 
be achieved in realistic experimental setups, as shown in 
the next section. 



C. Loop Entropy and Hook factors 

The general rates defined in Sec. |TT] include both the 
entropy loss factor and the hook factor. Both depend on 
the length q of the bubble, and for both their relative 
influence diminishes for increasing bubble lengths. In 
the Fokker-Planck description both factors arc omitted, 
which are the assumptions (ii) and (Hi). These assump- 
tions arc valid for long bubbles, which can be obtained 
by having long soft zones and keeping the temperature 
far above the melting temperature of the soft zones, i.e. 
Us > 1. 

Fig. [H] studies the effect of the loop exponent c on the 
coalescence time density 7r(t) for a fixed length of the 
barrier region, A^ = 15 and varying the lengths of the 
soft zones. Even for soft zones of length Nl = Nji = 30 
there is a substantial influence of the loop entropy fac- 
tor, making long soft zones a requirement in experimen- 
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Figure 6: Coalescence time density for varying soft zone 
lengths and fixed length of the barrier A'^ — 25. Included 
are plots for Ub — 0.98 (upper) and Ub — 1.1 (lower), i.e. for 
temperatures below and above the melting temperature of the 
barrier, respectively. Furthermore, c = fi = 0. 
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Figure 7: Coalescence time density for A'^ = 25, Nl ~ Nr 
20, Ub — 0.98 and different values of Us. Furthermore, c 
fj, = 0. 



tal realizations which should agree reasonably with the 
Fokker-Planck approach. 

The hook factor leads to a decrease in the rate constant 
k, k ^ kq~^^ where q is the length of a given bubble, so 
introducing the hook exponent /i > leads to a consider- 
able and bubble-length-dependent decrease of the tran- 
sition rates, see Fig. [101 However, for sufficiently large 
bubbles, the rate is roughly constant and introducing a 
renormalized rate constant k ~ kL^ can compensate for 
this effect. Here L is a characteristic bubble size. If the 
temperature is kept well above the melting temperature 
of the soft zones, and the length of the soft zones are 
much longer than the barrier, L is well approximated by 



the length of the soft zones, L 



N 



L/R- 



In Fig. [TO] we 



illustrate the effect of a renormalized rate constant (with 
L = 53 being the best fit) together with the standard 
rate coefficient k = 1. 

In conclusion, both the influence of the loop entropy 
factor and the hook exponent can be eliminated by keep- 
ing the length of the soft zones sufficiently long and using 
a renormalized value for the rate constant k. 



VIII. RELEVANCE FOR SINGLE MOLECULE 
EXPERIMENTS 



Relevant for the separation of statistical weights are 
the empirical relations Q 



rpAT 

m 

rpGC 



(355.55 + 7.95 ln[Na+])X (95a) 
(391.55 + 4.98 ln[Na+])X, (95b) 



which give the melting temperatures of GC and AT pairs 
in terms of the (intermediate) salt-concentration in the 
solvent obtained from melting experiments [65j . Note 
that the value of stems from an average over all 

possible combinations of AT and TA base pairs Q; if 
only TA/AT and AT/TA pairs arc interchangeably used 
the value of T^^ can be lowered further. The above re- 
lations can be translated into free-energy difference s by 
AG = AS'(T„-r), where AS" = -24.85cal/(mol K) 
Eqs. (|95ap . (j95bp contain contributions from both base 
stacking and hydrogen bonding, and is thus the melt- 
ing temperature suitable for our situation. It has been 
shown that the de pen dence on salt-concentration lies in 
the stacking term [66| and not as previously thought in 
the hydrogen bonding term [gJ ■ The stacking is a combi- 
nation of hydrophobic, electrostatic (screening of the neg- 
atively charged phosphate groups), and dispersive inter- 
actions but there is no app arent consensus on which term 
is the dominant one [Ba]- At high salt-concentrations 
~ 1 — 5M the temperature dependence levels off due to 
a decrease in the hydrophobic effect; with most water 
molecules tied up in the solvation of ions, the entropy 
decrease involved in base stacking is small [6l| . 

The melting temperature of AT bonds has a stronger 
dependence on salt-concentration, so we can increase the 
ratio uat/ugc by decreasing the salt-concentration. A 
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Figure 8: Single trajectories for the two zipper forks, based on the Gillespie (Monte Carlo) algorithm presented in App. [C] 
The values of ut and Us axe stated in the figures. The other parameters are Nl ~ Nr = 20, A'^ = 25, c — fi = 0. The dashed 
horizontal lines mark the boundaries between the barrier region and the soft zones, and the arrows indicate the ends of the 
sampled trajectories. 
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Figure 9; Coalescence time density n{t) for fixed barrier width 
= 15, a nd va rying the length of the soft zones , Nl = Nr = 
10 in Fig. |9(a)| and Nl = Nr ^ 30 in Fig. |9(b)[ respectively. 
The other parameters are U(, = 0.98, /i = 0. 
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Figure 10: Coalescence time density 7r(t) for non-vanishing 
hook exponent fi = 0.588, Nl ^ Nr = 50 and iV = 25 for 
the barrier case, mj, = 0.98. The results are shown with and 
without the renormalized rate constant k, where we have used 
k = 53"'^*^ ~ 10.3. The other parameters are Us = 10 and 
c = 0. 



further benefit is a lowering of the melting temperatures 
thus enabling experiments well below T^'^ ~ 100°C, 
which is the case when [Na^] = O.IM, the standard 
concentration in electrophoresis experiments. High 
temperatures have practical disadvantages such as 
formation of air bubbles and increased evaporation of 
solvent molecules (68j . 

Most relevant experiments on DNA have been con- 
ducted at ~ O.IM salt-concentration. Those specifically 
looking at the salt dependence of the melting temper- 
ature work in the range ^ 0.01 — IM. A conservative 
estimate of [Na"''] = O.OIM gives uqc ~ 1 and uat — 6 
at r ~ 95°C, which is sufhcient for the Fokker-Planck 
approximation to be valid. 

Concerning the possible length of a DNA construct it 
should be reasonable to work with segments up to 100 — 
200 bps. In a setup combining fluorescence correlation 
spectroscopy and fluorescence quenching, as introduced 
in Ref. [1^|, the DNA is free to diffuse around in the 
solution. In this case the limiting factor is the time it 
takes for the quencher to diffuse in and out of the confocal 
volume. 

Practically the experimental method of choice may be 
a dual optical tweezers setup in which the DNA con- 
struct, via some handles of double-stranded DNA, is con- 
nected to two beads held in place by the tweezers. While 
this allows to keep the DNA construct in place the force 
exerted on the chain is relatively small. However this 
would allow direct observation of the construct avoiding 
diffusional correction. The centre of the barrier region 
could be decorated with either a fluorophore-quencher 
pair, or markers such as quantum dots or small gold 
beads that can be visualized by a microscope. The influ- 
ence of the attached markers should decrease with longer 
barrier length. Having this setup in a flow cell the system 
could be triggered by flushing in a solution with either 
different temperature or salt concentration. This can be 
done relatively quickly [6^ . Once the two initial bubbles 
are thereby created it should be possible to measure the 
coalescence time calculated herein. Repeating the ex- 
periment would produce the distribution of coalescence 
times, from which important system parameters can be 
inferred. Decorating the barrier region with several, suf- 
ficiently small, markers would, in principle, allow one to 
measure the position of the coalescence. 



IX. FINAL CONCLUSIONS 

Single molecule techniques give us increasing insight 
into the behavior, equilibrium and dynamic, of biopoly- 
mers. Of particular and outstanding interest is DNA, due 
to its importance in biological contexts as well as its role 
as a model biopolymer. In order to extend our knowl- 
edge about the biological function of DNA it is crucial 
to quantify and understand the denaturation behavior 
of DNA at the single molecule level, its sequence depen- 
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dence and, ultimately, its relevance to genetic processes 
such as transcription initiation. A major question hereby 
concerns the dynamics of transient DNA denaturation 
bubbles. 

While first single molecule fluorescence correlation ex- 
periments have demonstrated the feasibility of monitor- 
ing the fluctuations of a single bubble, some questions 
remain about the model system and the explicit setup 
used in these experiments. In particular, the obtained 
time scales for base pair zipping and unzipping as well as 
the influence of the attached fluorophorc-quencher pair 
remain under debate. 

Here we suggest an alternative model for accessing 
DNA stability parameters and basepair (un-zipping) con- 
stants: in our model system DNA bubbles in two AT-rich 
regions are formed and separated by a more stable GC- 
rich barrier region. The coalescence behavior of the two 
DNA bubbles across the barrier region is then studied. 
We show that the stability parameters for bubble and 
barrier regions can indeed be chosen sufficiently differ- 
ent to allow preparation of the DNA construct in the 
proposed fashion, by the proper adjustment of tempera- 
ture and/or salt concentration. Once coalesced the newly 
created single bubble is stabilized against immediate re- 
closure of the barrier region both dynamically and due 
to the release of the boundary free energy corresponding 
to one cooperativity factor ctq. Appropriate fluorophore- 
quencher tagging of the barrier region basepairs should 
therefore allow for the direct observation of the bubble 
merging dynamics. 

Apart from the relevance of the investigated system 
for understanding the dynamics of DNA and its bio- 
logical function, the mathematical description presented 
here is of interest for its own sake as it corresponds to 
a previously unsolved case of two vicious random walk- 
ers in opposite linear potentials. We established the 
solution of this problem by solving a bivariate Fokker- 
Planck equation (continuum limit of the discrete master 
equation description) analytically. In a careful analysis 
we showed under what conditions the Fokkcr-Planck ap- 
proach is valid and what deviations one would expect for 
realistic systems. Furthermore, the analytic results were 
explained using qualitative arguments and corroborated 
using stochastic simulations. 
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Appendix A: CALCULATION OF p(x;t\xo) VIA 
LAPLACE TRANSFORM 



In this appendix we present a detailed calculation 
of the single-walker auxiliary density p{x;t\xo) satis- 
fying Eq. ((49|) together with the boundary conditions 
Eqs. (|37ap and (j37b[) . as well as the initial condition 
p{x;t — 0\xq) ~ S(x — xq). p{x;t\xo) solves the 
Schrodinger equation 



d ^, . . 
—p[x;t\xo} 



(92 

— -f 



p{x;t\xQ), (Al) 



that, after a Laplace transform and some rearrangement 
becomes 



(9^ 
ox'' 



p{x; z\xq) = -S{x - Xq), 



(A2) 



with k{z) = yz + f^ (wc skip the explicit z-dependence 
in the formulas from now on). 

Consider first the solution of equation 



dx' 



g{x,xo) = 0, 



(A3) 



for X < Xq with boundary condition (|37ap . and for x > Xq 
with boundary condition (j37bp . The solutions are 



g{x,xo) 



a 



X < Xq, 



X > Xq, 



(A4) 



with K= {k + f)/{k- f). 

The solution of Eq. (|A2[) can now be found by the 
ansatz 



p{x]z\xq) 



C (Ke^'' + e-*^'^) (e*=^« + Ke^'^e-'^^o) 



C (kc'^^" 



(A5) 



») (e''^ -t- ne^^e ''^) , x > xq, 



such that 

p{x] z|xo) - C {e^i^+^°) + Ke^'^e"'^'!^^^''!) 

(x+x,,)^ (A6) 



^k\x — X{)\ 



where C is determined from the jump condition by inte- 
grating Eq. ((X2|l 



— 1 = lim 

£^0 + 



xo+c r q2 



d 



/ 



^0+ ox 



p{x; z\xo)dx 
d 



xo+c dx 



p{x; z\xo) 



= 2fcC(l - K^e^'') . (A7) 
Here, it has been used that p{x; z\xo) is continuous. 
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Figure 11: Enumeration scheme for the numerical analysis: 
The two-dimensional grid points {xL,m) are replaced by a 
one- dimensional running variable s. See text for details. 



Appendix B: IMPLEMENTATION OF THE 
DISCRETE MASTER EQUATION 



To solve the eigenvalue equation (|86p by a numerical 
scheme, it is convenient to replace the two-dimensional 
grid points (x^,™) by a one-dimensional coordinate s 
counting all lattice points, compare with We choose 
the enumeration illustrated in figure [TT] From this figure 
we notice that m € and xl & [0, A/ — m]. An 

arbitrary s-point can be obtained from a specific {xl,iti) 
according to: 



s={m-l)M --^ '-^ -+XL. (Bl) 



From this relation we notice that the maximum s value 
is 



S = : 



c{s} = 



+ 1) 



1, 



(B2) 



i.e., the size of the relevant VF- matrix (see below) scales 
as M^/2. Expression (jBip allows us to change the trans- 
fer coefficients to the s-variable, t^^j^{xL,'ni) —y t^^j^{s), 
using the explicit expressions ([SO]) . (|8ip . and ([55)) for 
the transfer coefficients, together with the boundary con- 
ditions in equations ([TS]) and (|79p . From equation (|Bip 
and Fig.[Tl]we notice that [compare with Eq. ([86]) ] 



im+l 

\xL-l 
mi— 1 

\XL + 1 



s|™ + M - m, for xl > 1, 
s|™ - (M - m + 1), for m > 2, 
- (M - TO + 2), for TO > 2, 

s|^'^ +A/-TO+1, 

for XL < M - (to + 1) & to < M - 1, 



(B3) 

Eq. (|86p can then be written in matrix form as 

W{s, s')Qpis') = -VpQpis), (B4) 



where explicitly the matrix-elements are 
W{s, s + M ~ to) 



W{s,s- [Af-TO+1]) 
W{s,s-[M -m + 2]) 
W{s, s + M -m + 1) 



tj(.s + M - to), 
for s rtl a;L > 1, 
t^(.s- [M-to + 1]), 
for s rtl TO > 2 
t+{s- [A/-m + 2]), 
for s rtl TO > 2, 
<^(.s + A/-m + l), 
for s rtl XL < AI - (m + 1) 
& TO < A/ - 1, 

+t+{s)+t],is)), (B5) 



and the remaining matrix elements are equal to zero. 
We have introduced the notation s rtl with the meaning 
"s is to be taken for". The problem at hand is that 
of determining the eigenvalues and eigenvectors of the 
(5 + 1) X {S + l)-matrix W above. The coalescence time 
density is then calculated from Eqs. ([M]) and ([M|) . In 
terms of the running variable s, see Eq. (|Bip . and the 
VF-matrix defined in equation (jB5P the detailed balance 
conditions (fTOj) and (fTTj) become 



W{s,s')^{s') = Wis',s)3f{s). 
The orthogonality relation, Eq. ([87]) . becomes 
Qp{s)Qp'{s) _ 



(B6) 



(B7) 



Convenient checks of the numerical results then include: 
(i) The eigenvalues should be real and negative (so that 
rjp > 0); (ii) The eigenvectors should satisfy the orthonor- 
mality relation, Eq. (|B7[) . 



Appendix C: STOCHASTIC SIMULATION OF 
BUBBLE COALESCENCE 



In this Section we give a brief introduction to the 
stochastic simulation of DNA-breathing, for details we 
refer to Ref. [sl]. We apply the Gillespie algorithm in- 
troduced in 1976 as a stochastic approach to the study 
of chemical reactions [T^I • 

Following the schematic of Fig. [1] we simulate the dy- 
namics of the two zipping forks separating the two initial 
bubble domains from the barrier region. As each fork can 
either zip or unzip, the system is described by the four 
different rates, where fi G {+,—}, and v G {L,R}. 
Given these rates, we assume that the statistical weight 
for a given event, {fi,!^}, to occur in a time interval 
[t,t + 5t] is f^St. Then the idea of the Gillespie scheme is 
the following [t^ : The probability that nothing happens 
in the time interval [t,t + r], and that in the following 
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interval [t + T^t + t + dr] an event of type {fi, i/} occurs, 
is the so-called reaction probability density 



PiT,^i,ly)dT = Po{T)t'idT. 



(CI) 



To determine the probability Po{t) that no event hap- 
pens within [t, t + r], this interval is divided into K spans 
of duration e = t/K. The probability that no event oc- 
curs in the first subintcrval [t, t + e] is then 



ll[l-t':;e] = l-J2t^^e + 0{e'). 



(C2) 



IJL,V 



Treating the remaining intervals similarly produces an 
expression for Pq, 



Poir) 



-\ K 



l-Y^t^^e + Oie') 



(C3) 



Taking the limit K ^ oo and reinserting in Eq. (jCip . wc 
find the Poissonian law 



PiT,fx,iy)^t^cx_j>i-J20 



(C4) 



At some given instant of time, t, the system is in a cer- 
tain configuration. The update is performed as follows: 

(i) The rates are calculated according to the configu- 
ration. 

(ii) A set of random numbers (r, ^, h'), distributed accord- 
ing to P{t, ^, ly) in Eq. (jC4p . is drawn from a generator. 

(iii) The time is advanced according to t ^ t + t, and 
the configuration is updated according to the randomly 
chosen event /i, i'. 

The steps (i)-(iii) are repeated until a specified stop crite- 
rion is fulfilled, in our case the merging of the two initial 
bubbles. We record the stop time and the final configu- 
ration, and a new run is initiated using the same initial 
condition. 



Following Ref. [70| we briefly present how random 
numbers r and fi can be constructed using numbers 
drawn from a uniform distribution: Let PciT') be some 
continuous probability density function, e.g., Pc(T')dT 
is the probability for finding a r within the interval 
[r', r' + dr]. The associated probability distribution fimc- 
tion is then defined as 



Fciro) 



Pc{T')dT\ 



(C5) 



which is the probability of some r being less than tq. To 
get a random r according to Pc, given some random num- 
ber i? S [Oj 1] drawn from the uniform distribution, we 
have to invert Fc{t) = R. Using Pc{t) = y P{t, fi, i^) 
from Eq. (|C4[) . with r > and inverting the expression, 
we obtain 



In 



(C6) 



Similarly, we determine the appropriate random number 
for the direction of the "reaction" (zipping/unzipping of 
left /right zipper fork), following 



(C7) 



is the probability of having fj. < fiQ. Inversion given some 
random number R € [0, 1] drawn from the uniform dis- 
tribution is now requiring that Fd{fJ. — 1) < R < Fd{p-). 
Using Pditi) = j P{t, fj)dT the random event ^ is deter- 
mined by 



N 



(C8) 
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